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This paper derives an efficient procedure for using the three-dimensional (3D) vector radiative 
transfer equation (VRTE) to adjust atmosphere and surface properties and improve their fit 
with multi-angle/multi-pixel radiometric and polarimetric measurements of scattered sunlight. 
The proposed adjoint method uses the 3D VRTE to compute the measurement misfit function 
and the adjoint 3D VRTE to compute its gradient with respect to all unknown parameters. In 
the remote sensing problems of interest, the scalar-valued misfit function quantifies agreement 
with data as a function of atmosphere and surface properties, and its gradient guides the search 
through this parameter space. Remote sensing of the atmosphere and surface in a three- 
dimensional region may require thousands of unloiown parameters and millions of data points. 
Many approaches would require calls to the 3D VRTE solver in proportion to the number of 
unknown parameters or measurements. To avoid this issue of scale, we focus on computing the 
gradient of the misfit function as an alternative to the Jacobian of the measurement operator. 
The resulting adjoint method provides a way to adjust 3D atmosphere and surface properties 
with only two calls to the 3D VRTE solver for each spectral channel, regardless of the number of 
retrieval parameters, measurement view angles or pkels. This gives a procedure for adjusting 
atmosphere and surface parameters that will scale to the large problems of 3D remote sensing. 
For certain types of multi-angle/multi-pixel polarimetric measurements, this encourages the 
development of a new class of three-dimensional retrieval algorithms with more flexible 
parametrizations of spatial heterogeneity, less reliance on data screening procedures, and 
improved coverage in terms of the resolved physical processes in the Earth's atmosphere. 

© 2014 Elsevier Ltd. All rights reserved. 


1. Introduction 

More accurate and complete monitoring of cloud and 
aerosol properties is needed to reduce uncertainties in both 
the radiative forcing of climate and feedbacks between the 
radiative forcing and changes in global temperature that are 
the result of changes to clouds and their properties [1]. 
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http://dx.doi.Org/10.1016/j.jqsrt.2014.03.030 
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While multi-angle polarimetric measurements and plane 
parallel retrieval methods provide the capabilities necessary 
for regions that are horizontally homogeneous [2-4], the 
retrieval of aerosols in broken cloud fields and near cloud 
edges remains an open challenge which limits the coverage 
and accuracy of retrievals [5-7], Using the three-dimensional 
(3D) vector radiative transfer equation (VRTE) can address 
this issue by explicitly accounting for the spatial distribution 
of solar illumination, scattering material, and polarimetric 
measurements. In contrast to plane-parallel and spherical 
models for radiative transfer, the 3D VRTE places no default 
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Nomenclature 


Vectors are assigned to bold lowercase letters, 
matrices are assigned to bold uppercase let- 
ters, and integral operators are assigned to 
script uppercase letters. 

Remote sensing problem 

n, N index for unknowns and total number 
m, M index for measurements and total number 
L number of wavelengths 

a, a" vector of unknowns and elements 

y,y"' measurements vector and elements 

y(a),y"'(a) measurement model vector and elements 


0(a) 

Domain 

X 

h(x) 

Vh(x) 

V 

D 
dD 
D X 

r+ 

r_ 


misfit function 


position in space 

signed distance to boundary 

outward pointing normal 

direction unit vector 

unit sphere 

spatial domain 

spatial domain boundary 

internal set 

outgoing set 

incoming set 


Single scattering 

(7(x; a) volume extinction coefficient 
Z(x,v,v';a) volume scattering kernel 
R(x, V _ , V + ; a) surface reflection kernel 

Forward vector radiative transfer 


Pi 

9o 

Aqo 

n* 

ul 


polarization analyzer (internal adjoint source) 
polarization analyzer (outgoing adjoint source) 
internal source vector for misfit gradient 
outgoing source vector for misfit gradient 
adjoint of scattering operator 
adjoint of reflection operator 
adjoint solution operator 


Forward integral equations 


T"oo 

foi-ward streaming 
internal) 

operator (internal-to- 

r.o 

forward streaming 
internal) 

operator (incoming-to- 

ro+ 

forward streaming 
outgoing) 

operator (internal-to- 


forward streaming 
outgoing) 

operator (incoming-to- 

f 

internal solution to the forward integral 
equations 

g 

incoming solution 
equations 

to the forward integral 

Adjoint integral equations 


00 

adjoint streaming 
internal) 

operator (internal-to- 

-0 

adjoint streaming 
incoming) 

operator (internal-to- 

^ 0 + 

adjoint streaming 
internal) 

operator (outgoing-to- 


adjoint streaming 
incoming) 

operator (outgoing-to- 

p 

internal solution to the adjoint integral 
equations 

q 

outgoing solution 
equations 

to the adjoint integral 


u 

Stokes vector solution 

Subscripts 

/o 

internal light source (internal forward source) 



So 

solar illumination (incoming forward source) 

o 

indicates a source vector, i.e. right-hand side 


internal source for parameter derivatives 

0 

defined on or related to the internal set 


incoming source for parameter derivatives 

-1- 

defined on or related to the outgoing set 

Z 

scattering operator 

- 

defined on or related to the incoming set 

n 

reflection operator 



lAa 

forward solution operator 

Superscripts 

Adjoint vector radiative transfer 

T 

transpose of a vector or matrix 



❖ 

adjoint of an operator 

w 

adjoint Stokes vector solution 

\ " 

integration variables 


restrictions on the spatial variability of the atmosphere and 
surface [8], Work to extend coverage with 3D methods 
has shown promise for determining average cloud optical 
thickness [9] and cloud top height [10], However, as a side 
effect of the increased flexibility, the 3D VRTE leads to 
retrieval problems with many more unknown parameters 


and multi-pixel measurement constraints. A significant con- 
cern is therefore the extent to which a proposed algorithm 
scales “gracefully” to large problems. The objective of this 
work is to formulate an adjoint method for the 3D VRTE 
which maintains the scalability required for the application 
to atmospheric remote sensing problems. 


70 


W. Martin et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 144 (2014) 68-85 


Adjoint methods and other linearization procedures can 
reduce the number of radiative transfer simulations needed 
over the course of an iterative procedure for fitting data 
[11,12], During each iteration, the current estimate of the 
atmosphere and surface properties is adjusted to improve its 
fit with measurements. Making the right adjustment requi- 
res knowledge of how the fit will change when the adjust- 
ment is made. This in turn entails solving the 3D VRTE for 
each wavelength to evaluate the misfit between model and 
measurements and its gradient with respect to unknown 
parameters. It is worth noting that, brute-force numerical 
differencing could be used to compute the gradient with 
0(LN) radiative transfer computations, where N is the num- 
ber of parameters and L is the number of wavelengths. 
Adjoint methods provide an alternative route to computing 
this derivative, and analogous work in the field of medical 
imaging shows that the number of calls to the radiative 
transfer solver can be as low as two: one forward solve for 
evaluating the misfit with data and one adjoint solve 
for computing the gradient with respect to retrieval para- 
meters [13-15]. Although this result was specific to angularly 
averaged measurements in the frequency domain with wave 
length-independent parameters, it is an example of a scal- 
able adjoint method derived for an analogous scalar trans- 
port equation. We define a retrieval adjustment procedure as 
scalable if it can be applied to problems with arbitrarily many 
measurement constraints and unknown parameters without 
requiring additional calls to a radiative transfer simulation. 
A scalable method can require 0(L) radiative transfer solu- 
tions at each step, but not 0{LN) or 0(M), where M is the 
number of measurements (including all wavelengths, view 
angles and pixels). 

Previous work on adjoint/linearization methods for remote 
sensing provides this kind of improved efficiency for plane- 
parallel [16], spherical [17-19], and pseudo-spherical radiative 
transfer models [20-22]. Also, adjoint techniques have been 
used to approximate solutions to the 3D VRTE for atmospheric 
properties with small deviations from plane-parallel symme- 
try [23,24]. To our Imowledge, the use of adjoint methods to 
develop a scalable remote sensing methodology that relies on 
the 3D VRTE has not hitherto been considered. This topic is of 
current interest due to the advancement of computational 
techniques, which simulate 3D scalar and vector RT in the 
Earth's atmosphere [25-29], including the recently released 
vectorized Spherical Harmonic Discrete Ordinate Method 
code by Evans [30]. So far. the primary application of these 
codes has been synthetic studies which assess errors asso- 
ciated with plane-parallel retrievals [31,32]. In our view an 
equally important direction of research deals with how to 
incorporate simulations of the 3D VRTE directly into cloud 
and aerosol retrieval algorithms [9]; for example, using multi- 
pixel methods to improve single-pixel retrievals by accounting 
for adjacency effects, or perhaps, using a futuristic 3D para- 
metrization of clouds and aerosols to retrieve their spatial 
variability in complex scenes with broken cloud cover. The 
objective of this paper is to provide the necessary theoretical 
foundation for such endeavors, by extending adjoint methods 
to allow scalable computations of the misfit function and its 
gradient using codes that solve the 3D VRTE. 

To ensure that the adjoint method derived here meets 
the needs of the atmospheric remote sensing community. 


the theoretical description of measurements is consistent 
with ground-based, air-borne, and space-borne polarimeters, 
and the parameter-adjustment methods are similar to those 
used in operational retrieval algorithms [16,33]. The method 
focuses on minimizing a misfit function for passive measure- 
ments of scattered sunlight, but active measurements and 
measurements at other wavelengths may be included as 
prior constraints on spatial variability: using high spectral 
resolution LIDAR to constrain the aerosol scattering coeffi- 
cient [34] or microwave cloud tomography to constrain 
cloud-droplet volume concentration [35,36]. Moreover, we 
formulate the adjoint framework in a manner that is con- 
sistent with the complex microphysical parametrizations 
needed to model single-scattering properties in the Earth's 
atmosphere [37]. The procedure is outlined using the stan- 
dard integro-differential form of the 3D VRTE and derived 
using an equivalent integral formulation, written using 
concise operator notation for the processes of streaming, 
scattering and reflection. The integral formulation and asso- 
ciated operators are related to existing numerical solutions, 
and we describe how to extend such codes to solve adjoint 
radiative transfer by using the reciprocity principle to write 
the adjoint Stokes-vector solution in terms of a slight 
modification to the usual forward solution. 

Preliminary definitions are organized in Section 2 with 
the fundamental adjoint property asserted but left tempora- 
rily unproven. The use of adjoint methods in developing a 
scalable procedure for adjusting parameters as a part of a 
remote sensing methodology is described in Section 3. Then, 
the general framework for forward and adjoint 3D VRTE is 
derived in Section 4 and the fundamental adjoint property is 
proven in Theorem 1. Supplementary technical results are 
presented in Appendices A and B. 

2. Preliminaries 

This section introduces the theoretical framework of the 
forward and adjoint 3D VRTE as needed for large scale 3D 
remote sensing of the atmosphere and surface. We generalize 
adjoint methods to arbitrary boundaiy conditions, and this 
entails a mild reformulation of the forward 3D VRTE as a 
boundary value problem and proof of the fundamental adjoint 
property for the corresponding boundai-y value problem of 
adjoint 3D VRTE. The logical progression that we choose to 
follow is to define the forward and adjoint 3D VRTEs as inde- 
pendent boundary value problems. Then, we prove the funda- 
mental adjoint propeity which relates them. The first task is 
to define the domain. 

Definition 1 (Domain). The spatial region of interest is an 
open, connected, and bounded set D c with smooth 
boundary dD c — smooth to guarantee that the bounding 
surface has a continuous outward-pointing normal vector. 
The region is described implicitly by its signed-distance-to- 
boundary function: 

- inf llx-x'll forxG(DU4D) 

X' € dD 

inf lix-x'll for X G R^\(D U 4D). 

X' e dD 

The useful properties of h are that it is continuous on and 
differentiable near the boundaiy, with gradient Vh(x) equal 
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to the unit normal vector pointing out of the domain for each 
X e dD. The function value determines if a given point, x g 
is inside the region of interest. h{x) < 0; on its boundary, 
h{x) = 0; or outside. h(x) > 0. 

Taking the direction vector v to be always in the unit 
sphere, c IR^ ; we define the three regions making up 
the domain of the 3D VRTE: the internal set. 


D X = {(X, V): h(X) < 0}, 


(2) 

the outgoing set. 



r+ ={{x,v+ ): fi(x) = 0 and v+ 

■ Vh(x) > 0), 

(3) 

and the incoming set. 



r_ = {(x,v_): h(x) = 0 and v ■ 

Vh(x) < 0). 

(4) 


Inner products of Stokes-vector and source-vector func- 
tions are defined for each domain: the internal Inner 
product, 

")dxs2 = f [ dSv w(x, V) ■ u(x, V), (5) 

Jh{x) < 0 Js^ 

the outgoing inner product. 


for all square-integrable functions, /. The smallest such 
value. C. is called the operator norm, ||£||op, and it will be 
needed in Section 4,2 to state the constraints on scattering 
and reflection that guarantee solve-ability of the 3D VRTE. 
The adjoint of an operator is defined in the usual sense, as 
the operator. £*, which satisfies the adjoint property: 

{pX\f]) = {C*\plf}- (13) 

The adjoint, C*, gives the alternative rule for evaluating the 
inner product in Eq. (13), so that numerical procedures 
may use whichever side is more efficient. 

In summary, the three distinct subdomains for 3D vector 
radiative transfer are defined through the utility function, 
h(x)-. the interior set, D x s^: outgoing set, r+: and the 
incoming set. r_. Each subdomain has an inner product and 
a set of square-integrable functions. As a convention, we 
reserve operator notation for bounded linear operators to 
ensure similarity to matrix algebra. Lastly, the adjoint of a 
linear operator was defined. 

2.1. Forward and adjoint 3D VRTEs 


(q, u)r^ = f dSx [ dSv^ |v+ 

Jh(x) = 0 Jv+ Vh{x)>0 

■Vli(x)| q(x,v+)-u(x,v+), (6) 

and the incoming inner product, 

(W,g)r_ = [ dSx [ dSv_ |v_ 

Jh{x) = 0 Jv_-Vh(x)<0 

■ Vli(x)| w(x,v_) -g(x,v_). (7) 

These elementary inner products appear so often in pairs 
that is helpful to write the forward inner product as 


**IdxS^ 


= iPX)n^s^+{q,U)r 


( 8 ) 


DxS^©r+ 


and the adjoint inner product as 
m/Idxs^ \ if 

Ms,, 

' ' Dxs^ec- 


= ('^=/>DxS^ + (W,g)r 


(9) 


Each of the three inner products in Eqs. (5)-(7) defines a 
vector space of square-integrable functions, for example 
the internal source vectors, /, such that 

llfllLs^ = (fJ>DxS^<~. (10) 

or the incoming source vectors, g, such that 

llgllr_ =(g,g>r_ <oo. (11) 


Erom the operator point of view, the square-integrable 
functions in Eqs. (10) and (11) are vectors in a linear space, 
and linear operators will act in much the same way that 
matrices do. To guarantee this, we give symbolic repre- 
sentation only to linear operators that are bounded. By this 
convention, a linear operator with symbol. C, will act on a 
square-integrable vector. /, and return another square- 
integrable vector, £[f]. This follows from the definition of a 
bounded operator: the linear operator, L, is bounded if 
there exists a value. C, so that 

IK[f]ll<C|lf||, (12) 


The purpose of this paper is to formulate an efficient 
procedure for adjusting unknown parameters as a part 
of an abstract remote sensing problem. The atmosphere 
and surface properties are described by an unknown 
parameter vector, a = (a") for 0<n<N, from which 
physical single-scattering properties are derived: 
extinction, (T(x;a); scattering kernel. Z(x,v,v';a); and 
reflection kernel. R(x,v ,v+;a). Eor a given illumina- 
tion defined by incoming and internal light sources, the 
3D VRTE provides a solution vector u(x,v;a) to be used 
in modeling each polarimetric measurement as an 
inner product over a pair of polarimetric response 
functions, Pg and q™ : 

y'"(a) = {p™,u)o^s2-E(q^,u)r^, (14) 



with measurement vector y = (y'") for 0<m<M. This 
motivates the boundary value problem of the 3D VRTE 
which defines the Stokes vector solution, u. for incom- 
ing solar energy. 

Definition 2 (Forward 3D VRTE). for a fixed parameter 
vector, a, and the corresponding single-scattering proper- 
ties, a, Z, and R\ the forward solution u(x, v; a) is defined 
for square-integrable source vectors. /g and gg. as the 
unique solution to the integro-differential equations of the 
forward 3D VRTE: 

V ■ Vu-|-(tu-Z[u] =/g on D X S^, (16) 

u|^ -■7^[u|^^]=go onr_. (17) 

The integral operator for scattering is defined as 

Z[u](x,v) = ]- [ dSv-Z(x,v,v')-u(x,v'), (18) 

for (X, V) G D X S^, and the integral operator for reflection is 
defined as 
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K[ulrJ(x,v^) 



dSv^. |v+ ■ Vh(x)|R(x, v_,v+) ■ u(x,v+), 

•V/l(J£) > 0 

(19) 


for (x,v_)er„. Coupling with boundary conditions is 
imposed by spatial continuity along outgoing and incom- 
ing directions: 


ulrjx,v+) = limulg^g2(x-tv+,v+), 


( 20 ) 


u|;- (x,v_) = limu|( 3 ^ 52 (x-^tv_,v_), 


( 21 ) 


for (X, v_^) s r+ and (x,v^)er_. 

Assuming existence and uniqueness for the moment, we 
define the solution operator for the boundary value 
problem of the forward 3D VRTE, Ifal'}: 


{ 


«l 

«l 


DxS^ 
C + 



( 22 ) 


satisfy the adjoint property: 



(25) 


Notice that all solver operations now act on the fixed adjoint 
source vectors, ^Pq and Aqg , so that the change in fit can be 
computed for any unknown parameter by integration. This 
motivates the definition of the boundaiy value problem for 
the adjoint 3D VRTE. 

Definition 3 (Adjoint 3D VRTE). Eor fixed parameter, a, 
and single-scattering properties, a, Z, and R\ we define the 
adjoint solution, w, for square-integrable adjoint source 
vectors, Pg and q^, as the unique solution to the adjoint 
3D VRTE: 

-V ■ Vw-|-(7W-2*[w] =Pg on D X S^, (26) 


The forward solution operator, Ua, is a 2x2 matrix of 
integral operators which is parametrized by a and acts on 
internal and incoming source vectors, /g and gg, to give 
the Stokes vector solution on the internal and outgoing 
sets, and . An explicit formula for the forward 

solution operator is derived in Section 4.2.1 and given by 
Eq. (112). 

The solution operator, Ua, plays the role of a forward 
solver in the present discussion of an abstract remote 
sensing problem — one call to a forward solver is equiva- 
lent to evaluating the solution operator, Ua, for one pair of 
source vectors, /g and q^. Using this operator we write 
the model for polarimetric measurements in Eq. (15) as the 
inner product of detector response functions, p™ and q™ , 
with the Stokes vector solution for solar illumination: 


y"'{a) = 




DxS^mr+ 


(23) 


This expression allows the computation of all measure- 
ments with one call to the solution operator, Ua, followed 
by relatively inexpensive integrations over the polari- 
metric response function for each detector. Therefore, Eq. 
(23) is well suited to the task of evaluating many different 
measurements for fixed internal and incoming source 
vectors. 

However, in remote sensing applications where unlcnown 
atmosphere and surface parameters, a, are adjusted to fit data, 
the computation also requires variation of source functions for 
computing the components of the gradient of the misfit 
function, 30 1 da". We will see in Section 3 that the desired 
quantity has the following form: 


30 

da" 


APo 

Mg 



A/"g 

M"g 


^ DxS^®r+ 


(24) 


for fixed adjoint source vectors, Apg and Aqg, and forward 
source vectors, 4fo and Agg, corresponding to derivatives 
with respect to each parameter, for 0 < n < N. Eorward and 
adjoint source vectors appearing in Eq. (24) are defined 
explicitly by Eqs. (42), (43), (46), and (47). To avoid repeated 
calls to the solver, Ua, we seek the adjoint of the forward 
solution operator. This adjoint operator, (Ua)*, is defined to 


-7^*[wfy ] = qg onr+. (27) 

The adjoint-scattering operator is defined as 

2'*[w](x, V) = ^ [ dSv' Z(x, v', v)^ ■ w(x, v'), (28) 

for (x,v) gD X S^, and the adjoint-reflection operator as 

1t*[w\r_](x,v+) = ^ [ dSv_ |v_ ■ Vfi(x)|R(x,v_,v+)’' w(x,v_), 

(29) 

for (x,v+)gT+. Coupling with boundary conditions is 
imposed by spatial continuity along outgoing and incom- 
ing directions: 

wfy^(x,v+) = limw|o>,s2(x-tv+,v+), (30) 


wfy (x,v_) = limw|o,<s 2 (x-htv_,v_), (31) 

t\o 

for(x,v+)Gr+ and (x,v^)Gr_. 

Assuming existence and uniqueness for the moment, we 
conclude by defining the adjoint solution operator. Eor 
each parameter a, the adjoint solution operator, U*, maps 
adjoint source vectors, Pg and qg, to the adjoint Stokes 
vector on the internal and incoming sets: 



The explicit form of this 2x2 matrix of integral operators 
is derived in Section 4.2.2 and given by Eq. (119). 


Since the boundary value problems for the forward and 
adjoint 3D VRTE are defined independently, the adjoint 
property, (Ua)*=U*, requires proof, and this is done in 
Theorem 1 . In summary the forward 3D VRTE is stated in 
Definition 2 and can be used to evaluate M radiometric 
measurements with 0(L) calls to the forward solution 
operator, Ua. We asserted that the adjustment of atmo- 
sphere and surface properties would require evaluation of 
the left hand side of Eq. (25), and noted that this would 
require 0(LN) calls to Ua- In Section 3, we show how the 
left hand side of Eq. (25) arises naturally as the required 
quantity in an iterative search, and how evaluating with 
the adjoint alternative on the right-hand side of Eq. (25) 
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leads to a scalable procedure for adjusting unknown para- 
meters as a part of a remote sensing algorithm. That is, one 
which requires only 0(L) evaluations of the solution opera- 
tors at each step, making the number of calls independent 
of the size of the problem. 

3. Application to remote sensing 

In the context of remote sensing of the Earth's 3D atmo- 
sphere and surface, the adjoint method provides a means of 
adjusting arbitrarily many atmosphere and surface para- 
meters to improve their fit with arbitrarily many polarimetric 
measurements without changing the number of 3D VRTE 
simulations needed. To provide a concrete example of this, 
consider the task of retrieving cloud, aerosol, and surface 
properties for Yellowstone National Park which has a surface 
area of ten thousand l<m^. For measurements, suppose we 
have access to a hundred satellite images of the park — taken 
with a single-spectral channel, from different perspectives, 
and with 1 km resolution. These data provide one million 
constraints. M = 1 x 10®. Suppose also that a discretization of 
the atmosphere and surface is constmcted with a total of one- 
thousand volume and surface elements, and that the volume 
single-scattering and surface reflection properties are repre- 
sented by an average of ten parameters per discrete element. 
These rough assumptions would result in the use of ten thou- 
sand parameters to describe the cloud, aerosol, and surface 
properties. N = 1 x 10^. 

At each step in the retrieval algorithm we must adjust 
these ten thousand parameters to decrease the collective 
misfit with one million measurements. We note that if one 
were to linearize the measurement operator for this problem 
then the Jacobian matrix, consisting of elements dy^/da", 
would have ten billion entries. One of the key strategies of the 
method outlined here is to avoid the computation and storage 
of the Jacobian matrix, working instead with the misfit 
function and its gradient. Since the misfit function is scalar 
valued, its gradient in this case has only ten thousand 
elements. Even in this extreme example, the adjoint method 
described here provides a procedure for adjusting parameters 
to improve the collective fit with all data using only two calls 
to the 3D VRTE solver. For multi-wavelength data the required 
number of calls is 0{L). 

Although this example describes the scalability of the 
adjoint method using a futuristic application involving a full 
3D reconstmction of cloud and aerosol properties, it is worth 
noting that this fits within a hierarchy of methods that start 
with retrievals assuming a plane-parallel atmosphere and 
with each pixel being an independent column [16,38,39]. This 
approach has been extended by Dubovik et al. [33] to include 
statistical modeling of the co-variation of atmospheric and 
surface properties in different pixels within the framework of 
a multi-pixel optimal estimation scheme. A natural addition 
to their usage of a multi-pixel prior probability distribution 
would be the usage of a multi-pbcel measurement operator, in 
which the 3D VRTE couples the radiative effects of nearby 
columns and allows, in the context of clear sl<y observations, 
for the proper account of adjacency effects. In this context, the 
adjoint method would provide a means of adjusting plane- 
parallel retrievals to correct for 3D, or adjacency effects. 
Moreover, the scalability result implies that the number of 


calls to the 3D VRTE solver is independent of the number of 
columns or pixels so that adjustments can be made to many 
pixels at once. 

The remainder of Section 3 will summarize the methodol- 
ogy which makes this scalability possible. Qualitatively the 
adjoint method accomplishes this by associating the residual 
misfit between model and measurements with a single source 
distribution for the adjoint 3D VRTE. The residual for each 
individual image pixel is defined as the difference between 
model and observation, and is specific to the location of the 
instrument, the field of view of the pixel, and the sensitivity of 
the polarization analyzer. The weighted sum of these localized 
and directed residuals over all image pixels gives a single 
distribution of adjoint sources, and the adjoint 3D VRTE is 
solved to back-propagate this residual through all orders of 
multiple scattering. Then, simple integrals can be evaluated to 
determine the change in fit for all possible adjustments to the 
unknown parameters. This alternative way of thinldng pro- 
vides a rule for computing the misfit gradient with the desired 
scalability. The subsequent use of the misfit gradient in 
numerical optimization routines is discussed in Section 3.3, 
where each iterative adjustment to cloud, aerosol, and surface 
properties is written as a solution to an JVx N system of linear 
equations. 

3.1. Model for polarimetric measurements 

The chosen setting for these results is the Earth's 3D 
atmospheric shell bounded between the Earth-atmosphere 
interface and an arbitrarily large radius out to space. The 
results extend to any smooth connected sub-region of inter- 
est, provided that reasonable horizontal boundary conditions 
are imposed. In the context of the radiative transfer model of 
light propagation, incoming solar radiation is scattered in the 
atmosphere and reflected by the surface, causing measur- 
able radiative effects that vary with location, direction, and 
polarization. 

To setup an abstract remote sensing problem, lety = (j/™) 
for 0 < m < M be a vector of single-wavelength multi-angle/ 
multi-pixel polarimetiic data taken by a ground, air. or space 
borne instrument. Let a = (a") for 0 < n < JV be a vector of N 
unlcnown parameters which define a three-dimensional dis- 
tribution of cloud, aerosol, and surface properties. In practice 
there will be constraints on the parameter vector, a, to gua- 
rantee a reasonable physical interpretation (e.g. non-negative 
particle concentrations). Enforcing such constraints by finitely 
many linear equalities and convex inequalities is ideal, to give 
a numerically convenient description of the convex set of all 
possible states of the atmosphere and surface. We now 
describe in three steps, how the definitions of Section 2 lead 
to a useful model for polarimetric measurements of atmo- 
spheric radiation as they depend on cloud, aerosol, and 
surface properties. 

First, the values of the volume-extinction coefficient, 
a(x;a); the volume-scattering matrix, Z(x,v,v';a); and the 
surface- reflection matrix. R(x,v ,v+; a), must be written 
explicitly as smooth functions of the vector of parameters, a. 
They must be smooth to guarantee the existence of deriva- 
tives, dajda", dZjda", and dRIda": and also to guarantee 
that the integral operators defined in Eqs. (46) and (47) 
will return square-integrable functions. Furthermore, for all 
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feasible values of the parameter vector, the single-scattering 
properties must satisfy solve-ability criteria given in Section 4 
by Eq. (109). The parametrization of single-scattering proper- 
ties incorporates both spatial and micro-physical variability. 
Surface reflection at the Earth-atmosphere boundary is 
characterized by several spatially dependent parameters. 
Volume-extinction and scattering properties are modeled as 
a linear combination of contributions from molecular scatter- 
ing and various modes of airborne-particle. For each mode 
there are parameters that define loading, size distribution, 
shape, and complex refractive index. For the present discus- 
sion. we assume that there is a well-defined functional 
relationship between parameters and single-scattering prop- 
erties. that the functions are smooth with respect to para- 
meters. and that the single-scattering properties satisfy the 
solve-ability criteria. 

Second, the 3D VRTE in Definition 2 is used to solve for 
the multiple scattering of incoming solar radiation and 
determine the Stokes vector solution. u{x,v,a). For each 
feasible parameter vector, a, the method requires a solver. 
Ua, that acts on forward source vectors for solar illumina- 
tion. /g and gg. and returns the Stokes vector solution 
(restricted to the internal and outgoing sets): 


{ 


«l 

u| 


DxS^ 

r+ 



(33) 


The restrictions. u|£,.,g 2 and . of the full Stokes vector 
solution, u, provide all the information that is necessary to 
model multi-angle polarimetric measurements. 

The third step involves expressing the measurable quan- 
tities which correspond to elements of the data vector, y, as 
inner products of the solution with detector response func- 
tions. Internal measurements are computed as the inner 
product of the internal Stokes vector, with a polar- 

ization analyzer defined on the internal set. Pg(x, v): 

y'"(a)=(p™,u)( 3 ^s 2 forO<m<Mi. (34) 

Outgoing measurements are computed as the inner product 
of the outgoing Stokes vector. u|,- . with a polarization 
analyzer defined on the outgoing set. qg(x,v+): 

y’^(a)={q1,u\r^)r^ for Mi <m<M. (35) 

Given the clear apertures of typical Earth observing instru- 
ments we note that the polarization analyzers will be effec- 
tively Dirac-delta dishibutions in the location variable with 
angular integrations being determined by the field of view of 
the given sensor or pixel. While Dirac distributions are not 
square-integrable functions, they may be approximated as 
such to within discretization error. Aircraft measurements 
taken inside the domain result in a weight that is localized to 
a point in space: 

p1(x,v)ocS(x"'-x). (36) 


For ground based instruments, e.g. AERONET [40]. a natural 
route to computing measurements is to localize the position 
of the instrument to a point on the boundaiy: 


q™(x+,v+)cc 


5(x™-x+) 
iv+ ■ vfi(x™)r 


(37) 


Satellite measurements may be taken at a great distance away 
from the domain, so in this case we suggest projecting the 


data to the outgoing boundary. This results in a weighting 
distribution that is singular in direction: 


q™(x+,v+)cc 



|v+ ■ Vfi(x+)| 


(38) 


Again, this singular distribution can be integrated over an 
actual instmment field of view in practice, with a scaling by 
the reciprocal. |v+ ■ V/r(x+)| to counteract the weight that 
appears in the inner product Internal measurements 

require no such scaling. In this way. the formalism used is 
shown to be consistent with common types of radiometric 
and polarimetric measurements taken of the atmosphere. 

To summarize the process of modeling polarimetric mea- 
surements the three steps are as follows: (1) Single-scattering 
properties are written as smooth functions of N parameters 
that describe cloud, aerosol, and surface properties. (2) The 
Stokes vector solving the 3D VRTE is computed as a model for 
the spatially and directionally dependent field of radiative 
energy in the atmosphere. (3) Each individual polarimetric 
measurement is represented as an inner product of the 
internal or outgoing solution with a polarimetric response 
function, defined on the same domain. This procedure pro- 
vides the theoretical connection between observations and 
the retrieval target of atmospheric composition. 


3.2. Data misfit and gradient calculation 


In the abstract remote sensing problem, we aim to use 
multi-angle polarimetric data stored in the M-dimensional 
vector, y, to adjust the 3D atmosphere and surface para- 
meters stored in the N-dimensional vector, a, and reduce 
the measurement residual, y-y(a), to within measure- 
ment error. Using the instrument's measurement error 
covariance matrix. Se, the misfit of model and data is 
quantified by the value of the misfit function: 

= ^(y -y(a))'^ -Sf^ ■ {y -y(a)) . (39) 


To improve the fit we seek to adjust unknown parameters. 
a, to decrease the value of the misfit function. The steepest 
decrease in cf> is obtained locally in the direction opposite 
to the gradient. V0. The nth element of this JV-dimensional 
vector is given by the following formula: 


d0(a) 

da" 


(y-y(a)Y 


dy(a) 
da" ' 


(40) 


By differentiating Eqs. (34) and (35) with respect to 
parameter a" and collecting terms into integration kernels 
for the internal and outgoing data, we write Eq. (40) as a 
sum of inner products: 


d0 

da" 



dU \ 
do"/ 



(41) 


The adjoint source vectors for differentiating the misfit 
function. ^Pq and Aqg. are obtained by summing over all 
detector response functions: 

Apg(x,v;a) =2: E (y"''-y"''(a))(Sf^)m'mPo(x,v), 

0 <m' < MO < m < M] 

(42) 
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and 

Aqg(x,v+;a) 

0 < m' < MMi <m<M 

(43) 


However, the alternative rule for computing the gradi- 
ent with the adjoint 3D VRTE is analytically equivalent to 
Eq. (49), but requires only one additional call to a 3D VRTE 
simulation, independent of how many parameters are 
used. The adjoint rule for differentiating the misfit func- 
tion is as follows: 


These adjoint source vectors, Apg and Aq^,, may be 
visualized as collections of many “search lights” emanating 
from all measurement pixels at once, with the intensity of 
each search light equal to a weighted sum of the measure- 
ment residuals. 

The next key step is to evaluate the derivative of the Stokes 
vector, du/(5a". Although this could be accomplished numeri- 
cally by finite-difference methods, a better way is to solve the 
3D VRTE with modified volume-source and incoming-Stokes 
vectors. To see how this is possible, we differentiate Eqs. 
(16) and (17) with respect to parameter a" and obtain the 
following 3D VRTE for du/da": 


30 



APo 

Mo 


A/o 


^ DxS^®r_ 


(50) 


Using the fundamental adjoint property. (Uaf =U*, which 
is proven in Section 4.3 as Theorem 1, along with Eq. (32); 
the gradient can be written in terms of the adjoint Stokes 
vector, w, solving the adjoint 3D VRTE with adjoint source 
vectors, APo and Aqgj : 



(51) 


^dU 

dU 


3a" 

3a" 

dU 


dU 


r -Ti 


3a" 


da" 


dU 

3^ 


= A/" 


= Ag"o, 


(44) 

(45) 


with right-hand sides equal to fonward source vectors, A/^ 
and Ag((, . The internal source vector, A/g (x, v) for (x, v) s D x 
S^, accounts for the change in extinction and scattering: 


\fl(x,v,a) 


(X; a)u(x, V; a) 

3an^ , ) y ^ ) 

^ f dZ 

dSv'^x,vy-a)u(x,v';a). 

4;t7§2 dQ^ 


(46) 


And, the incoming source vector, Agg (x, v_ ) for (x, v_ ) s r_, 
accounts for the change in reflection: 

Ag^(x,v_;a) =^[ dSv+ |v+ 

V/i(x) > 0 

■ Vfi(x)| ^(x, v_ , v+ ; a) • u(x, v+ ; a). (47) 


Comparing with Definition 2, we see that the left-hand side of 
the 3D VRTE for the gradient of the Stokes vector solution, 
3u/3a", is identical to the left-hand side of the 3D VRTE of the 
solution itself, u. Therefore the same existence and unique- 
ness results are applicable. If the forward source vectors for 
computing parameter derivatives are square-integrable then 
the solution operator returns the derivative of the Stokes 
vector with respect to the parameter, a": 


( 3U \ 



\ 3a" 

> =Ua- 

A/"g 1 

— 1 
[ 3a" 

0 = 1 


Substituting Eq. (48) into Eq. (41 ) we can express the equation 
for the gradient of the misfit function as 


-^=(W,Arg)o,s2 + (W,Ag”g)r_. (52) 

The significance of Eqs. (51) and (52) is that the adjoint 
solution, w, is independent of which a" appears in the diffe- 
rentiation. Therefore any component of the gradient can be 
evaluated with the same adjoint solution, w, and without 
further calls to a radiative transfer solver. The procedure 
involves the comparatively inexpensive operations of weigh- 
ting the adjoint solution with the source terms from Eqs. 
(46) and (47) and integrating over the internal and incoming 
sets. The computation time of these integrations is (at worst) 
comparable to a single order-of-scattering computation, and 
could be much faster if one is careful to use a sparse basis for 
spatial and directional variability. 

Using the adjoint method to compute the gradient of the 
misfit function shifts all the multiple-scattering computations 
to the residual distribution which depends only on the current 
atmospheric state and the misfit between observations and 
the polarized radiance that is generated by the current 
atmospheric state. The number of solutions to the 3D VRTE 
required in evaluating Eqs. (51) and (52) is 0(L) and inde- 
pendent of the number of parameters, N. Therefore, mles for 
adjusting cloud, aerosol, and surface parameters based on the 
adjoint calculation of the misfit gradient are scalable to 
retrieval problems with many measurements and unknown 
parameters. 


3.3. Iterative parameter adjustment 

To discuss how the adjoint computations of the gradient of 
the misfit function can be incorporated into a scalable 
retrieval algorithm, we define a regularized misfit function, 

0Teg- 

0reg(a) = 0(a) + 0pnor(a)- (53) 


3a" ~\\ Ago 



(49) 


for each a" with 0 < n < N. As written in Eq. (49). computing 
all elements of the gradient requires N solutions to a 3D VRTE 
solver at each step of a multi-step iterative procedure. 


Prior information is introduced by ^pnor, to give a new maxi- 
mum-likelihood estimation problem which is less sensitive to 
measurement noise and to provide a means of imposing 
additional measurement constraints, for example those from 
a coordinated LIDAR instrument [41,42]. The retrieval starts 
with an initial guess, Oq, and makes additive adjustments, bi^. 
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SO that the updated parameter, = a^ + b);, converges to a 
minimizer of the regularized misfit function, cf>reg. 

A common starting place for many optimization meth- 
ods is the second order Taylor approximation for ^reg, 
about any feasible state, a: 

^reg(tt-t-b) 

= a>reg(a)-|-va)reg(a) • b + ^b^ ■ VV0reg(a) ■ b. (54) 


which approximates the Hessian of the measurement 
misfit function using the gradient of the misfit function, 
V0, at previous iterations. These gradient-based methods 
can take advantage of the scalability of the adjoint rule in 
Eq. (50). The Hessian of the misfit function is approxi- 
mated using previous parameter estimates, a^, and pre- 
vious gradients. V0(aif). The approximate Hessian. is 
updated and improved upon at each step: 


With sufficient prior information the Hessian matrix, 
VV0reg(fl). can be made positive definite. In this case the 
minimum value of the Taylor approximation can be found 
by solving Newton's equation for step, b: 

(VV<P(a)-rVV<Pprior(a)) • b= -(V0(a)-rV0prior(a)). (55) 


A common scenario in setting up Newton's equations is 
that the Hessian of the prior function is easy to compute 
and the Hessian of the measurement misfit function is 
prohibitively expensive. Quasi-Newton methods substitute 
an approximate Hessian. 

In atmospheric remote sensing applications, methods 
such as Levenberg-Marquardt use a linearized measure- 
ment model to approximate the Hessian in terms of the 
Jacobian matrix [33,39,41 ]. The Jacobian matrix. 


J(aQ 


dy(ak) 

da 


(56) 


is computed at every step, to approximate the Hessian as 
follows: 


WW0{a)^J{af -m. 


(57) 


V V0(afc) Hkiao, . . ., Ofc, V0(ao), V0(a^)). (60) 

The rule usually guarantees symmetry and non-negative 
definiteness, and this is the case for the Broyden-Fletcher- 
Goldfarb-and-Shanno method, which is used in the med- 
ical imaging applications [13-15]. With a positive definite 
prior Hessian matrix VVcPprior, the approximate Newton's 
equation uniquely defines a parameter adjustment, b^, 
with the following linear system: 

(f/fc-h VV0pi.ior(flk)) b= -(V<Z>(ak)-hV0prior(a)<)). (61) 

Provided that the problem is well scaled and that the 
approximate Hessian, is chosen appropriately, the step 
will result in an improved set of cloud, aerosol and surface 
properties via the updated parameter =a);-i-bj(. More- 
over, by using the gradient, V^, to set up the local problem, 
the method can leverage the scalability of the adjoint com- 
putation to adjust atmosphere and surface properties with 
only 0(1) calls to a 3D VRTE solver at each step. 

3.4. Pseudo-forward problem 


The Jacobian matrix contains the derivatives of measure- 
ments with respect to unknown parameters and is, in 
general, dense with 0(MN) elements. Using the funda- 
mental adjoint property, the inner product for each ele- 
ment may be written in either of the two forms: one using 
the forward solver. 


dy'^ 

~da^ 


fAn 

ro ’"“1 


DxS^ ®T4 


and the other using the adjoint solver, 

A/; 


da" 


= (U\ 


P'e 


Ag"o 


DxS‘'ec_ 


(58) 


(59) 


This section describes how to solve the adjoint 3D VRTE 
using a computer simulation for solving the forward 3D 
VRTE. To do this we must define two actions, a and Q, 
which transform vectors according to the following rules. 
Action by a changes the sign of the direction argument: 

af(X,V)=f{X,-V). (62) 

Action by Q, flips the orientation of circular polarization: 


Qf(x,v) = 


- 

■/l(X,V)' 


/q(x,v) 


fu(x,v) 

_ 

fv(x,v) 


(63) 


In Eq. (58) the forward solution operator. Ua, must be 
evaluated for each different parameter and wavelength and 
this results in 0(LN) computations. If the number of para- 
meters is small or the computational solver provides Green's 
functions, then Eq. (58) can be quite efficient, see for example 
[19,22[, but codes which currently solve the 3D VRTE do not 
compute Green's functions. In the alternative mle given by Eq. 
(59) the adjoint solution operator, W*, must be evaluated for 
each measurement and this results in 0(M) computations. 
However, it is worth mentioning that alternative approaches 
use the single-scattering approximation to formulate an 
approximate, sparse Jacobian matrix. This idea has shown 
promise for retrieving the volume-scattering coefficient, a, 
using the scalar 3D radiative transfer equation with isotropic 
and wealdy scattering media [43]. 

To handle the large amount of multiple scattering in 
clouds, we considered another quasi-Newton method 


These actions are their own inverses. Qf = oi^ = Identity. 

For scattering media which obeys the principles of 
mirror-symmetry and reciprocity, the kernel for the 
single-scattering operator satisfies the rule. 

Z^(x, v', V) = QZ(x, - V, - v')Q, (64) 

and the kernel for the reflection operator satisfies the rule, 

R^(x,v',v) = QR(x, -V, -v')Q. (65) 

The adjoint scattering operator can be written in terms of 
the forward scattering operator, 

Z*[W] = a<lZ[aQwl (66) 

and the adjoint reflection operator can be written in terms 
of the forward reflection operator, 

7^*[W] = aQn[aQW\. 


(67) 
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By plugging Eqs. (66) and (67) into the adjoint 3D VRTE 
defined by Eqs. (26) and (27) and acting on both sides with 
aQ, it is easy to verify that the transformed adjoint Stokes 
vector, aQw, solves the forward 3D VRTE with pseudo- 
forward source vectors, «QPo and «Qqo. Therefore, the 
adjoint solution operator can be evaluated using the 
forward solution operator as follows: 


u 


* 

a 


{ 


Po 


90 


a QUa 


«QPo 1 

“0.90 J ■ 


(68) 


This means that a computer code that solves the forward 
problem can be used to solve the adjoint problem. Pro- 
vided that it is sufficiently general to accept the trans- 
formed source vectors, «QP0 and “090. it will output as a 
solution the transformed adjoint Stokes vector, aQw. The 
only additional difficulties in solving for the adjoint solu- 
tion, w, arise in preparing the right-hand side and inter- 
preting the solution. 


4. The fundamental adjoint property 

This section defines mathematical tools for proving the 
fundamental adjoint property with the proof given in 
Section 4.3. The first objective is to define a family of 
streaming operators. These will enable the formulation of 
integral equations which are equivalent to the integro- 
differential equations of 3D VRT, as given by Definitions 2 
and 3 in Section 2. The integral equations are presented in 
Sections 4.2.1 and 4.2.2, along with series expansions for the 
solution operators, Ua and U% Lastly, in Section 4.3 we state 
and prove the fundamental adjoint property as Theorem 1. 
This theorem shows that the adjoint 3D VRTE given in 
Definition 3 is well defined, that (Ua)* = U%, and justifies the 
use of this property in deriving a scalable procedure for 
adjusting 3D atmosphere and surface parameters. 


4.1. The streaming operators 

Streaming refers to propagation of radiative information 
along special line segments called chords.' Chords are defined 
to be the open-ended line segments in D whose endpoints lie 
on the boundary, dD. The streaming operators propagate 
source vectors along these chords: forward streaming opera- 
tors propagate sources in the positive direction, v, and adjoint 
streaming operators propagate sources in the negative direc- 
tion, -V. In contrast to previous derivations of the ancillaiy 
integral equation for 3D radiative transfer, for example that by 
Davis and Knyazikhin in Chapter 3 of [8], we split the strea- 
ming process into four distinct linear operations. This splitting 
facilitates the treatment of general boundary conditions in 
both the forward and adjoint 3D VRTE. Moreover, the 
streaming operators for the adjoint 3D VRTE are. actually, 
the adjoint operators of the forward streaming operators. 
A brief summary of the splitting of streaming operators will 
suffice for readers wishing to move ahead to Section 4.2 and 
the definition of the integral equations for 3D VRTE. 

Each streaming operator acts on, and returns, a func- 
tion which is defined on one of the three subdomains of 


’ Also called characteristics. 


vector radiative transfer. Eor instance, the internal stream- 
ing operator, T oo, acts on an internal source vector, /, 
defined on the internal set. D x S^. and it returns a Stokes 
vector restricted to the internal set, Boundary 

conditions are managed by the other operators. The 
incoming-to-internal streaming operator, T_o, acts on an 
incoming source vector, g, defined on the incoming set. r_, 
and it returns a Stokes vector on the internal set, 

The other two forward streaming operators are named 
according to their behavior in a similar way. The internal- 
to-outgoing streaming operator, Tq+, acts on an internal 
source vector, /, and returns an outgoing Stokes vector, 
: and the incoming-to-outgoing streaming operator, 
T_ + , acts on an incoming source vector, g, and returns an 
outgoing Stokes vector, . 

The adjoints of these streaming operators act on adjoint 
source vectors and return adjoint Stokes vectors. For 
instance, the adjoint of the internal streaming operator, 
Tqq, acts on an internal adjoint source vector, p, defined on 
the internal set, D x S^, and it returns an adjoint Stokes 
vector restricted to the internal set, w|p^g 2 . Flowever, the 
domains of input and output functions are reversed: the 
adjoint of the incoming-to-internal streaming operator, 
T*_o, acts on an internal-adjoint source vector, p, and 
returns an incoming-adjoint Stokes vector, w|,- ; the 
adjoint of the internal-to-outgoing streaming operator, 
acts on an outgoing-adjoint source vector, q, and 
returns an internal-adjoint Stokes vector, wl^^gi; and the 
adjoint of the incoming-to-outgoing streaming operator, 
T*__^, acts on an outgoing adjoint source vector, q, and 
returns an incoming adjoint Stokes vector, w|^ . 

Adjoint streaming operators are defined to satisfy the 


following adjoint properties: 


iP,T 00[f])DxS^ = ('^OoM J)dxS^> 

(69) 

(9,^o+[n>r, =<r*^[q],/>o,s2. 

(70) 

<P.^-0[g]>DxS^ = <^-o[P].g>C-, 

(71) 

(q,T.^\g])r,={T*_^[q],g)r_. 

(72) 


The remainder of Section 4.1 is devoted to parametrizing 
chords for the purpose of defining explicit rules for 
evaluating each of the forward and adjoint streaming 
operators. Forward streaming operators are defined in 
Eqs. (89)-(92), and their adjoints in Eqs. (93)-(96). In 
Appendix B we prove the adjoint properties of the stream- 
ing operators. 

4.1.1. Chords and boundary points 

Because streaming operators propagate information 
along chords, they are most easily defined with a chord 
parametrization. The signed distance-to-boundary func- 
tion, h(x), is quite useful for this purpose and enables 
treatment of non-convex atmospheric regions and surface 
topography. Using this function we define the unique 
chord for every internal, incoming, and outgoing point. 

For internal points, (x.v)eD x S^. we define the chord 
parameters as follows: 

t = xv. 


(73) 
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x^=x-tv. (74) 

Extreme values of chord parameter, t, are found by looking 
along the directions, v and -v, to the nearest boundary 
points: 

t_ = max{t_ G R: t_ < t and +t_ v) = 0}, (75) 

t+ = min{t+ G IR: t+ > t and h(x^ v) = 0}. (76) 

For outgoing points, (x, v+) g r+, the chord parameters 
are defined as follows: 

t+=x-v+, (77) 

x^ =x-t+v+. (78) 

The opposite extreme is found by looking along, -v+, to 
the nearest boundary point: 

t_ = max{t_ G R: t_ < t+ and h(x^ +t_ v+) = 0}. (79) 

For incoming points, (x, v^) g r_, the chord parameters 
are defined as follows: 

t_=xv_, (80) 

x^=x-t_v_. (81) 


measure zero in the integrals of interest and can therefore 
be neglected. We now focus on defining the streaming 
operators, noting that certain points that correspond to 
chords that are tangent to the boundary may require 
alternate definitions. 


4.1.2. Rules for evaluating streaming operators 

Streaming operators are defined by changing the argu- 
ment of evaluation from standard representation to chord 
representation, e.g. from (x,v) to (x^ -i-tv, v). This is done 
to isolate the direction, v, along which source vectors are 
integrated. The four forward streaming operators are 
defined as follows: 


T oo[fl(x, V) = r oo[f](x-^ + tv, V), 


= dt' exp^-^ dt" ct(x-^ - i-t"v)j/(x-'- -i-t'v, V) 


(89) 


To+ [f](x+ , v+ ) = To+ [f](x-L -h t+ V, V), 



df' 



dt" a(X^ -|-t"V) /(X^ +t'v,V) 


(90) 


The opposite extreme is found by looking along, v_ , to the 


nearest boundary point: 

t+ = min{t+ G R: t+ > t_ and h(x-^ -|-t+ v„) = 0}. (82) 

In each case, the position-direction pair can be asso- 
ciated with the unique non-empty chord, 

x'(t';x^, V) = x^ -l-t'v fort_<t'<t+, (83) 

through the interior, with h(x'(t')) < 0. The endpoints of 

the chord correspond to parameters, t_ and t+, 

x_ =x-^-l-t_VG(3D, (84) 

x+ =x-^-l-t+VG(3D, (85) 

and are on the boundary, since h(x_) = /r(x+) = 0. This 


defines the chord representation for all points in the 
domain of 3D VRTE and provides a useful alternative 
representation of internal, outgoing, and incoming points: 

(x,v) = (x-^-l-tv,v), (86) 

(x,v+) = (x^-ht+v+,v+), (87) 

(x,v„) = (x^ -l-t_v_,v_). (88) 

Looking toward future work implementing numerical 

methods, it is worth noting that this procedure can find 
chords through complex geometries by determining the 
zeros of a real-valued, single-variable function: h(x'f)). 

There is an important caveat. Although the chord is 
uniquely defined for each case, the endpoints, x_ and x+, 
do not always correspond to elements in the incoming or 
outgoing sets. The reason for this is that some chords will be 
tangent to the boundaiy at one or both end points. This 
poses a challenge to defining boundary-streaming operators, 
because there is not necessarily an incoming or outgoing 
point that corresponds to a location at which we desire to 
know the value of the streaming operator. However, the set 
of such points related to boundary-tangent chords will have 


^-o[g](x,v) = r_otel](x-^ +tv,v), 

rt 


= exp - 


dt" (tIx-*- -i-t"v) ]g(x-^ -i-t_ v,v), 


(91) 


r _ + [g](x+ , v+ ) = r _ + [g](x^ t+ V, V), 


= exp - 


(-/;■ 


dt" ct(x-'- -i-t"v) )g(x-'- -i-t_v, V). 


(92) 


The adjoints of these operators are given by the rules: 
V) = r*o[p](x-L -t tv, V), 


j: 


dt' 


exp( - dt" a(x^ +t"v) )p(x^ -i-t'v,v) 


(93) 


2'o+[q](x,v) = rS + [q](x^ +tv,v). 


= exp - 


(-r 


dt" cr(X-'- -|-t"V) )q(x^ -l-t + v, V), 


(94) 


2'*_olP](X- , v_ ) = r*_o[P](x-" + 1_ V, V), 


/; 


dt' 


exp( - dt" (j{x^ +t"v) )p(x^ -i-t'v,v) 


(95) 


T*_ ^ [q](x_ , ) = r* + [q](x^ -F t_ V, v), 

= exp^- ^ dt" (t(x-^ - l-t"v)^q(x^ -l-t+v, V). (96) 

The streaming operators in Eqs. (89), (90), (93), and (95) act 
on internal source vectors,/ orp, and require integration over 
some or all of the chords associated with the point of 
evaluation. Alternatively, the streaming operators in Eqs. 
(91), (92), (94), and (96) act on incoming forward source 
vectors, g. or outgoing adjoint source vectors, q. These involve 
only scaling by an attenuation factor. Note that while the 
operator Tq+ in Eq. (90) integrates over a chord from t_ to 
t+. the corresponding adjoint-streaming operator ^0+ in 
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Eq. (94) has no such integral. The integral over the chord is matrix-operator presentation is preferred over the use of 
subsumed in the inner product in Eq. (70). indices which are less easily readable. 


4A.3. Properties of streaming operators 

Streaming operators are defined to help solve the 3D 
VRTE and its adjoint by transformation to an equivalent 
system of integral equations. Under the action of the 
advective derivative, the forward streaming operators 
behave as follows: 

(v-V + <7)[Tooin]=f, (97) 

(V- V-hi7)[r_o[g]] = 0. (98) 

Similar properties hold for the adjoint streaming opera- 
tors: 


4.2.1. Forward integral equations 

We now present the integral formulation of the 3D 
VRTE and the so-called successive order of scattering 
series expansion for its solution. As described in the 
context of scalar radiative transfer by [25], the Stokes 
vector, u, can be written in terms of the solution vectors, 
/ and g, of the forward integral equations: 

“lDxS^=7'oo[f]+r_o[g], (104) 

U\r,=To+m + 'T- + \gl (105) 

U|r_=g- (106) 


i-V-V + a)[Tlo\p]]=P, 

(-V- V + i7)[rS^.[q]] = 0. 


(99) The vectors, / and g, are called solutions to the forward 
integral equations because they solve the integral formu- 
(100) lation of the forward 3D VRTE: 


These are shown in Appendix A as Theorem 2, and they 
provide the connection between integral and integro- 
differential forms of the 3D VRTE. 

4.2. Integral equations for 3D VRTE 

The integral operators for scattering and reflection. Z 
and TZ, act on Stokes vectors and return source vectors. The 
integral operators for streaming. Tqo,Tq+,T _o, and T_ + , 
act on source vectors and return Stokes vectors. In an 
approximate-numerical setting the Stokes vectors and 
source vectors will be represented by a finite-dimen- 
sional vector of basis function coefficients. Eurthermore, 
the integral operators will be approximated by linear- 
matrix transformations acting on these coefficient vectors. 
With the discrete analogue of matrix algebra in mind, we 
use a matrix operator notation to keep track of integral 
operations. As with the solution operators, Ua and U% curly 
brackets are used. {::} and {: }. The objective is to organize 
the linear integral operations according to the familiar 
notation of matrix- vector products from linear algebra: 

f .EToo ET.o'lf/'l _f ZToom + Zr.olg] 1 

\nTo+ “|7^ro+[f]+7^r_ + ter] J- 

(101) 

Normal array operations (such as associativity) behave as 
expected, for instance, the combined operations of 
streaming and then scattering can be written in either of 
the following two ways: 

f 2iToo ZT_o 1 f Z j f Too T_o 1 
[nTo+ 72r_+| = { 72 }|to+ T_+|- ^ ^ 

Empty spaces are assumed to represent a null operator. 
The analogy with matrices extends to the calculation of 
adjoints: 



(107) 


This differs from the ancillary integral equations (for 
diffuse radiation) as they are written in Chapter 3 of [8], 
in that Eq. (107) treats the internal source vector, /g, and 
incoming source vector, gg. as separate entities. This 
allows us to include direct radiation in the solution, 
causing the source vectors, /g and gg. to be identical to 
those of Definition 2. In Appendix A, Theorem 3, we show 
that the set of Eqs. (104)-(107) provide a solution to the 
3D VRTE that satisfies the integro-differential formulation 
in Definition 2. 

Eq. (107) is in standard form for a Fredholm integral 
equation of the second kind. To derive the series expansion 
for the solution operator, we use a fixed point iteration to 
obtain the following expression: 

f/1 ^ f ZToo ZT_o ^ I ZTao ZT_a 

\gj \KTo+ 7Zr_+J \gj TlT.+ ]\g^y 

(108) 

noting that powers, indicate repeated application of 
the integral operator. The solve-ability condition on the 3D 
VRTE must provide that the first term on the right-hand 
side of Eq. (108) will decay to zero. This occurs when a, Z, 
and R are such that the combined operations of streaming 
and scattering/reflection give an operator with norm less 
than unity: 


J TToo 
\l2To+ 


ZT_o 
72T_ + 



(109) 


If this condition is satisfied, we let K->oo in Eq. (108) to 
obtain the successive order of scattering series solution: 


f/1 ^ - f 2^00 Tr_oj"f/gj 

\sj 7^T_+J|ggJ■ 


( 110 ) 


f Too 
\To+ 



/t-H* 

0 



(103) 


where the adjoint of the 2x2 matrix-operator is defined 
using the joint inner products from Definition 1. This 


One may provide a more formal justification, as in [44], 
using the completeness of the space of square-integrable 
functions, as defined in Section 2. 

The series in Eq. (110) converges as a square-integrable 
function and provides the solution vectors, / and g, of the 
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foi-ward integral equations and by streaming them accord- 
ing to Eqs. (104) and (105), they provide the Stokes vector 
solution to the integro-differential form of 3D VRTE: 



We can thus define the solution operator for the forward 
3D VRTE, introduced in Definition 2: 



OO 

X I 

k = 0 


r_o 

T- + 

ZTqo 

7^ro+ 7^r_+j|goj■ 


( 112 ) 


4.2.2. Adjoint integral equations 

The adjoint 3D VRTE has a completely analogous 
integral formulation to that of the forward model. Using 
adjoint streaming operations, we write the adjoint Stokes 
vector, w, in terms of the solution vectors, p and q, of the 
adjoint integral equations: 

'VlDxs^=^oolP]+'^S+[qL (113) 

wir_=r*_o[p]+r*_ + m (114) 


W|r+ =q 


(115) 


The vectors, p and q, are called solutions to the adjoint 
integral equations because they solve the integral formu- 
lation of the adjoint 3D VRTE: 


Jpir^^rSo ^*^0+ _/Po 1 
jqj TZ*T*^j\qj jqQj- 


(116) 


The proof that the adjoint Stokes vector, w, given by these 
equations satisfies the adjoint 3D VRTE is given in Appendix 
A as Theorem 4. 

The vectors, p and q, can be expressed as the infinite 
successive-order of scattering series: 


fpl _ - I 


(117) 


This expression provides the solution vectors, p and q, of 
the adjoint integral equations and by streaming them 
according to Eqs. (113) and (114), they provide the Stokes 
vector solution to the integro-differential form of the 
adjoint 3D VRTE: 


{ 


w| 

w| 


DxS^ 




(118) 


The solution operator for the adjoint 3D VRTE, introduced 
in Definition 3, can therefore be written as follows: 



^ r z*T*o 1 '' r Pg 1 


(119) 


4.3. Fundamental adjoint property for the 3D VRTE 


Eor the solution operator of the adjoint problem in 
Definition 3, we wrote, “U*,” anticipating that it would be 
the adjoint of the solution operator for the forward VRTE, 
which we denoted by “(Ua)*". The notation hints that the 
equation, (Uaf =U% holds; that the adjoint of the solution 
operator for the forward 3D VRTE is the solution operator 
for the adjoint 3D VRTE. In fact, this relation justifies the 
name adjoint 3D VRTE, and must be proven through 
verification of the fundamental adjoint property: 



( 120 ) 

An equivalent form of this equation is written in terms of 
the solutions to the integral equations of 3D VRT, and it is 
this statement that we prove. 


Theorem 1 {Fundamental adjoint property). For any feasible 
parameter, a, and the corresponding single-scattering prop- 
erties, a, Z, and R, satisfying the solve-ability constraint in Eq. 
(109): the following fundamental adjoint property holds: 



where the vectors,/ and g, solve the forward integral equations 
with square-integrable source vector,/ ^ and ; and the vec- 
tors, p and q, solve the adjoint integral equations with square- 
integrable adjoint source vectors, Pg and q^ . 


Proof. We begin by substituting the adjoint integral 
equation for Pg and q^, given by Eq. (116), into the left 
hand side of Eq. (121) and proceed through the following 
steps: 



(122) 



r-+ 



DxS^®r+ 




00 

n + 1 

/T-* 

-0 

/' 

/ no 

n+ ' 

\no 

7'!.+ 

ZT 00 

ZT_ 

TZTo + 

7^r_ 




^ DxS^ ®r_ 



DxS ec- 


(123) 
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(124) 


These three steps are justified as follows: Eq. (122) is 
obtained by substitution of the integral equations of 
adjoint 3D VRT; Eq. (123) is obtained by using elementary 
adjoint properties of streaming, scattering and reflection 
operators; and Eq. (124) is obtained by substitution of the 
integral equations of forward 3D VRT. The second step can 
be verified by expanding to a sum of elementary inner 
products (that is (-.Odxs^- and applying 

elementary adjoint properties, and collapsing back into the 
matrix notation. □ 


5. Conclusion 

Adjoint methods can enable the use of 3D VRTE simula- 
tions for adjusting 3D atmospheric properties to fit multi- 
angle, multi-pixel polarimetric measurements of the Earth's 
atmosphere. This is shown by focusing on computing the 
misfit function and its gradient, and doing so with only 
two calls to a 3D VRTE solver for each wavelength. Scalable 
methods such as the adjoint method presented here will 
allow the role of the 3D VRTE to transition from a test bed for 
verifying plane-parallel retrievals to the core engine for perf- 
orming large-scale retrievals of atmospheric properties for 
scenes with strongly heterogeneous cloud cover. The primary 
benefit is that the 3D spatial dependencies of the sampling 
volume for remote sensing measurements will be explicitly 
modeled. The lack of default assumptions on cloud horizontal 
variability will allow for a more flexible parametrization of 
cloud stmcture and a more realistic model for measurements 
of broken cloud fields and the regions near cloud edges. As a 
near-term application of the adjoint method, Section 3 
discussed the use of a multi-pixel measurement operator to 
correct plane-parallel retrievals that have errors caused by 3D 
effects, including adjacency effects. Other applications are 
more ambitious and will require future research into how to 
parametrize cloud and aerosol properties in 3D. The imple- 
mentation of the adjoint method using existing codes is a 
topic of ongoing research, and we welcome collaborations in 
this effort. The adjoint method provides a way to adjust 
atmosphere and surface parameters that will scale to large 
problems — a foundation for a new class of retrieval algo- 
rithm. which uses a three-dimensional parametrization of 
atmosphere and surface properties to simultaneously recon- 
stmct both spatial and microphysical variability. 
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Appendix A. Equivalence of integral and differential 
equations of 3D VRT 

Integral and differential forms of 3D VRT were used 
interchangeably in the main text of the paper. This appendix 
Justifies such a use by showing that both the forward and 
adjoint system of integral equations provide a solution to the 
corresponding differential 3D VRTE. The first theorem states 
and proves properties from Eqs. (97) to (100). Then these 
properties will be used in the two theorems that follow: one 
for forward equivalence and another for adjoint equivalence. 

Theorem 2 {Streaming properties). The advective derivative 
acts on streaming operator output functions according to the 


following rules: 

(V- V-Ei7)[r_oteG] = 0, (A.1) 

{V-V + a)[Too\f]]=f, (A.2) 

(-V- V-Ef7)[r*+[q]] = 0, (A.3) 

(_v. V + i7)[r*o[p]]=p. (A.4) 


Proof. We begin with Eq. (A.l): 

V- vr_o[g](x,v) 

= ^exp^- ^ dr (7(x-^ -i-t"v)^g(x-^ -i-t_v,v), 

= -(7(x)r_o[g](x,v). 

Next we show Eq. (A.2) using an extension of the Leibniz 
rule to non-constant limits of integration [45]: 

V- VToo[f](x,v) 

= dt'exp^-^ dr (7(x^ -l-rv)^/(x-^ -l-t'v,v), 

= exp(0)/(x^ -I- tv, V) 

- (t(x^ -I- tv) dt' exp (^~ 

+ t"V)y(X^ +t'V,V), =/(X,V)-ff(X)Too[f](X, V). 
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For Eq. (A.3) we have the following: 

-V- vr*+[q](x,v) 

= -^exp^- ^ dt" cr(x-^ -i-t"v)^q(x^ -i-t+v,v), 
= -^7(x)r*^[q](x,v). 

Lastly, we show Eq. (A.4): 

-V- vrSo[p](x,v) 
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+ t"v)^p(x-^+t'v,v), 

= exp(0)p(x-^ + tv, V) 

- (7(x ^ + tv) dt' exp dt" a{x^ 

+ t"v)"jp{x^+t'v,v), =p(x, v)-ff(x)r*o[p](x,v). 

The four properties are thus verified. □ 

Theorem 3 (Forward 3D VRTE equivalence). The Stokes 
vector u given by the integral formulation o/3D VRT in Eqs. 

(104) -(107) solves the integro-differential 3D VRTE from 
Definition 2. 

Proof. First, we note that the functions u\r^ and u\r_ agree 
with the limits of u|(j^g 2 along lines approaching the bound- 
ary. The boundary conditions are verified by substituting Eqs. 

(105) and (106) into the left hand side of Eq. (17) and applying 
Eq. (107) to show the equality with the right-hand side: 

=g-7^[ro+[/^+r_ + [g]], 

=g-7^ro+[f]-F7^r_ + [g], 

= Sg- 

The solution on internal points is verified by substituting Eq. 
(104) into the left hand side of Eq. (16) and applying Theorem 
2 and Eq. (107) to show the equality with the right-hand side: 

V ■ VU + <tU-Z[U] 

= (V-S/ + a)[Toolf] + r_o\g]] 

-2[Too[n + r^o[g]], 

=f-ZToolf] + Z'r-o\gl 
=/o 

Thus, the Stokes vector u constructed from solutions / 
and g of the integral equations solves the differential 3D 
VRTE. □ 

Theorem 4 (Adjoint 3D VRTE equivalence). The adjoint 
Stokes vector, w, given by the integral formulation of the 
adjoint 3D VRTE in Eqs. (113)-(116) solves the integro- 
differential form of the adjoint 3D VRTE from Definition 3. 

Proof. First, we note that functions w|r_ and agree 
with the limits of w|p^g 2 along lines approaching the 
boundary. The boundary conditions are verified by sub- 
stituting Eqs. (114) and (115) into the left hand side of Eq. 
(27) and applying Eq. (116) to show the equality with the 
expected right-hand side: 

W\r^ -■7^*[W|;- ] 

= q-n*[T%m+r*_ + m, 

= q-n*T*_o\p]+n*T*_ ^ [q ] , 

= qo 

The solution on internal points is verified by substituting 
Eq. (113) into the left hand side of Eq. (26) and applying 
Theorem 2 and Eq. (116) to show the equality with the 
right-hand side: 

- V ■ VW + aW — Z*[W) 


= (-v- v+i7)[r*o[p]+rS+[q]] 

= p-z*r*olP]+z*r*^[q], =pg. 

Thus, the adjoint Stokes vector w constructed from the 
solutions p and q of the integral equations solves the 
differential 3D VRTE. □ 

Appendix B. Elementary adjoint property results 

Proving the adjoint properties for streaming operators 
will require us to equate certain multi-variable integrals, and 
this is facilitated by changing co-ordinates to integrate along 
chords. After summarizing these coordinate transformations, 
we will prove the four elementary adjoint properties for 
streaming operators stated in Eqs. (69)-(72). Eollowing this, 
we will show the adjoint properties of scattering and 
reflection operators. 

The natural basis to use for streaming integration proofs 
depends on the direction v s S^, and without loss of gen- 
erality we take the standard basis for v in terms of angles S 
and cp. This gives the following orthonormal basis for M^: 

V = [ sin 9 cos cp, sin 9 sin cp, cos 
5 = [cos 9 cos cp, cos 9 sin cp, — sini>]^, 
y» = [- sin y), cos cp, 0]^. 

The coordinate transformation for the internal set D x is 
given by the following rules: 

t = v X, 
y^ = 9 - X, 
y^=cpx, 

9 = 9, 
cp = cp. 

The rules change coordinates so that the chord parameter, t, 
controls the projection of x along v, while variables and 
y^ control location in the 2-dimensional plane orthogonal to 
V. We note that the perpendicular component used in the 
body of the paper, x^, is related to y’ and y^: 

x^ = 9y^ -h<py^. 

The associated surface element dy’ dy^ is written more 
compactly as dSx± . Eor evaluating integrals over the internal 
set, there is no change in weight associated with the chord 
transformation and 

dSv dSx± dt = dSv dx’ dx^ dx^. (B.l) 

For position direction pairs on the outgoing and incoming 
sets, the change of weight is exactly the cosine of the angle 
of incidence. For coordinate transformation from the per- 
pendicular component of a chord, x^, to the outgoing 
position, x+, we write 

dSv dS*± = dS*^ dSv^ |v+ ■ V/i(x+)|, (B.2) 

and for coordinate transformations to the incoming position, 
X _ , we write 

dSv dSx± =dS*_ dSv_ |v_ ■ V/i(x_)|. (B.3) 

These transformations are helpful in proving the four 
elementary adjoint properties for streaming operators. 
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Theorem 5 {Internal streaming adjoint property). The internal streaming operators T oo and T^q, defined by Eqs. (89) and (93), 
satisjy the adjoint property: 


(p,roo[n>D.s^ = (rSo[p]J>, 


DxS-^ 


(B.4) 


Proof. Beginning with the left hand side of Eq. (B.4), the equality is demonstrated in the following steps: 

[ dVx [ dSvp{x,vf ■Too\f](x,v)= [ dSv [ dS*x 2 f dt [ dt' p(x^ +tv,vfi- 
Jd Js^ Js^ Jx^(D) 4t_ 


exp - 


/ dsJ 


/'•dx( 

f^dt 

Js^ Jx^lD) 

(f _ ,f+ ) •• 

lt_ \q 

'r 


V)' 


expf - [ (7(x^+rv)dt"^p(x-L+tv,v)^ ■f(x^+t'v,v),= [ dS^ [ dS*x [ dt' +t'v, 

\ Jr / / Js^ Jx^iD) 

■f(x^+t'v,v), = f dVx [ dSvTlo\p](x,vf -fix.v}. 

Jd Js^ 

This is the right-hand side of Eq. (B.4). We note two key steps. Eirst, we rewrote the integral over space as an integral over 
chords: 


[dV,^[ dS,x 2 f* 

Jd Jx-^(D) Jt- 


dt. 


HD) (t-.t+) 

Second, we exchanged the order of integration of df and dt': 

^f+ ^t + 

L "''L 

Then, changing back to the original coordinates completed the proof. □ 

Theorem 6 {Internal source to outgoing Stokes vector streaming adjoint property). The streaming operator Tq+ defined by Eq. 
(90) has adjoint defined by Eq. (94): 


(q,^o+[f]>r, =(rS+[q],/>D,s^ 


(B.5) 


Proof. Beginning with the definition of the left hand side of Eq. (B.5) we proceed through the following steps: 
f dSx^ [ dSv+ |v+Vh(x+)|q(x+,v+)^-To+[f](x+,v+) = / dSv / dS*x +f+v. v)^ 

JdD Jv+-Vfl>0 Js^ Jx^lD) t+ 


dtexpl - 


a{x^ +t"v)dt"'^f(x^ +tv,v), = J^^dSv I ^^^dSxc'Z / 




-t"v)dr jq(x-L -rt+v,v)j f(x^ +tv,v), = J^dVxJ^^dSy [q](x, v)^ f(x. 


V). 


These steps show the equality. The key was writing each outgoing position x+ as the endpoint x+ = x^-i-t+vofa chord through 
the domain. Recognizing the definition of and changing to standard coordinates complete the proof □ 

Theorem 7 {Incoming to internal Stokes vector streaming adjoint property). The streaming operators T_g and T’Lg defined by 
Eqs. (91) and (95) are adjoint: 


(p,r _o[g])DxS^ = {T*_o\P],g)r 


(B.6) 


Proof. We begin with the definition of the left hand side of Eq. (B.6) and show the equality with the following steps: 

[ dVx [ dSvPix.vf ■T_o\g](x,v) = [ dSv/ dSxc'Z [ dt p(x^ +tv,vf ■ 

Jd Js^ Js^ Jx^iD) t- Jt- 

exp^- ^ cr(x-^ -l-rv)dr^g(x-^ -l-t_v,v), = dt 

expT- [ ct(x-l - i-t"v)dr^p(x-L -l-fv,v)^ -glx-L -l-t_v,v), = [ dS»_ [ 

\ Jt- ) ) JdD Jv_-Vh<0 

•g(X-,V_). □ 


dSv_ |v_ ■ V/r(x^)|r*_g[p](x_,v_)^ 
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Theorems (Incoming to outgoing Stokes vector streaming adjoint property). The streaming operators T - + andT*_^ defined 
by Eqs. (92) and (96) are adjoint: 

(q,T.^\g])r,={r*_^[qlg)r_. (B.7) 

Proof. We begin with the definition of the left hand side of Eq. (B.7) and show the equality with the following steps: 

f dSx / dSvjv+ ■ Vh(x+)|q(x+,v+)^ ■r„+[g](x+,v+), = / dSv / dS^c Y,q(x^ +t+v,v)'^- 

■IdD Jv+Vh>0 JS^ .Ih^(D) t + 

exp^-^ a(x^ +t"v)dt"^g(x^ +t^v,v), = J ^dSvJ dS^c Jj^exp^- o-(x-^ +t"v)df''^ 
q(x^ +t+v,v)') g(x^ +t_v,v), = [ dSx_ [ dSvJv_- 

/ JdD JV_Vh<0 

Vh(x„)|r*_ _^[q](x_,v_)^-g(x_,v_). □ 

Theorem 9 (Scattering operator adjoint property). Scattering operations Z and Z* are adjoint with respect to the internal inner 
product: 

(P, ^[«])dxs^ = “>DxS^ ■ (B-8) 

Proof. The left hand side of Eq. (B.8) is shown to equal the right-hand side by interchanging the order of integration over 

S": 

(p, Z[u])o>,s2 = [dVx [ dSv p(x, v)^ -j*- / dSr Z(x, v, v') ■ u(x, V), = [ dVx [ dSy- [ dSy^pix, vf ■ Z(x, v, v') 

JD Js^ Jo Js^ Js^ 

•u(x,v'), =^dV*^^dSv. (^^^^dSvZ(x,v,v'/-p(x,v)^ -u(x,v'), = (Z*[p],u)j^^g2. □ 


Theorem TO (Reflection operator adjoint property). The reflection operations TZ and TZ* are adjoint: 

(W|^ , 7^[Ur+ ]>r_ = (7^*[W|r ], Ur+ >r. 

Proof. The left hand side of Eq. (B.9) is equated with the right-hand side by interchanging the order of integration: 

(W|/- = [ dSx_ [ dSv_ |v_ ■ Vfi|w(x_,v_)’^ 

JdD 7v_V/l<0 


u. 


dSv, 


v+ ■ Vh 


R(x_,v_, v+) u(x_ 


,v+), = [ dS*_ [ 

JdD Jv+ 


dSv , |v+ ■ Vh| 


hi 


dSv_ |v_ ■ Vhl R(x_,v_,v+)^ ■ w(x_,v_) ■ u(x_,v+), = ],u| 


□ 


(B.9) 
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